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ABSTRACT 



We use high-resolution numerical simulations to study the physical properties 
of subhalos when they merge into their host halos. An improved algorithm 
is used to identify the subhalos. We then examine their spatial and velocity 
distributions in spherical and triaxial halo models. We find that the accre- 
tion of satellites preferentially occurs along the major axis and perpendicular 
to the spin axis of the host halo. Furthermore, the massive subhalos show 
a stronger preference to be accreted along the major axis of the host halo 
than the low-mass ones. Approximate fitting formulae are provided for the 
physical properties of subhalos. Combined with analytical and semi-analytic 
techniques, these empirical formulae provide a useful basis for studying the 
subsequent evolution of subhalos and satellite galaxies in their hosts. Future 
studies should however account for satellites that may not be undergoing the 
first infall in their evolution. 

Key words: large-scale structure of Universe - cosmology: theory - dark 
matter 



1 INTRODUCTION 



Cold Dark Matter (CDM) dominated cosmological models have received strong support from 
a wide range of observations (e.g., Spergel et al. 2003). In CDM models, the structures of 
dark matter are formed through the merger and accretion of smaller structures (White & 
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Rees 1978). Dark matter halos are high-density structures that are a few hundred times of 
the mean background density and are in approximate dynamical equilibrium according to 
the virial theorem. Galaxies are expected to condense in dark matter halos due to dissipative 
cooling processes. The observable properties of the galaxies are thus strongly influenced by 
the hierarchical growth of their host halos. During the hierarchical assembly of the halos, 
smaller halos merged into bigger ones are found to be long-lived even though their outer 
part may be stripped by the tidal interactions (e.g., Tormen 1997; Klypin et al. 1999; Moore 
et al. 1999). Some of these subhalos are believed to be the hosts of the satellite galaxies 
in galactic-sized halos (such as the Milky Way) and of the galaxies in clusters of galaxies, 
therefore it is essential to understand the evolution and dynamics of these subhalos in galaxy 
formation. Indeed, it has been demonstrated that it is important to resolve the merger of 
subhalos in order to accurately predict the cluster luminosity function of galaxies (Springel 
et al. 2001) and the luminosity function of field galaxies (Kang et al. 2004). The presence of 
the subhalos was used to interpret the anomalous flux ratios of lensed quasar images (e.g., 
Mao & Schneider 1998; Kochanek & Dalai 2004), though it is still unclear if the amount 
of subhalos predicted by the CDM models is in quantitative agreement with the lensing 
observations (Mao 2004; Mao et al. 2004; Zentner et al. 2005a). In addition, the subhalos 
have important effects on the heating of galactic disks (Benson et al. 2004) and possible 
7-ray emission from annihilation of dark matter particles (Taylor & Silk 2003; Stoehr et al. 
2003; see Bertone et al. 2004 for a review.) 

Because of the importance of the subhalos in cosmological studies, there have been many 
high-resolution N-body studies of subhalos since Moore et al. (1999) who discovered that 
subhalos can survive within bigger halos for a significant period of time. The mass function 
of subhalos was shown to be proportional to ~ (M sub /Mhost) _1 ' 8 (e.g., Klypin et al. 1999; 
Moore et al. 1999; Ghigna et al. 2000), and the radial distribution of subhalos above a given 
mass is much flatter than that of the dark matter in the host halo (e.g., Gao et al. 2004a; 
Diemand et al. 2004; Mao et al. 2004). The velocity dispersion of the subhalos is larger than 
that of the underlying dark matter (positive velocity bias), except in the innermost region 
where the velocity anti-bias is found for the subhalos (Ghigna et al. 1998; Okamoto & Habe 
1999; Colin et al. 2000; Klypin et al. 1999; Springel et al. 2001; Diemand et al. 2004). In their 
recent work, Gao et al. (2004a) showed that most present-day subhalos identified in their 
simulations have been accreted into their host halos very recently. This indicates that the 
subhalos accreted earlier, even if they were massive enough to host bright galaxies, may have 
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been disrupted in the current generation of high-resolution simulations (Gao et al. 2004b) 
due to physical (such as tidal stripping) and numerical effects (e.g., limited resolution). It 
is therefore very demanding resolution-wise to correctly model the dynamics, mergers and 
evolution of subhalos and satellites in the densest regions (Kang et al. 2004), and thus an 
analytical theory would be preferred to follow the evolution of subhalos within their hosts. 

Indeed, parallel to the numerical studies, significant progress has been made in under- 
standing the evolution and the distribution of subhalos with analytical models (e.g., van den 
Bosch et al. 1999; Klypin et al. 1999; Taylor & Babul 2002, 2004; Sheth 2003; Zentner & 
Bullock 2003; Ougri & Lee 2004; van den Bosch, Tormen & Giocoli 2005). The two principal 
physical processes, dynamical friction and tidal stripping, which determine the evolution 
of subhalos, are relatively well understood. However, as noted by Benson (2005), the ini- 
tial conditions (velocities and positions) of the infalling subhalos adopted in these studies 
are not from the prediction of hierarchical CDM models. Although some of these studies 
predicted subhalos with physical properties in reasonable agreement with those found in 
simulations, their success depends in part on the fine-tuning of the model parameters, such 
as the Coulomb logarithm function in the dynamical friction formula (e.g., Oguri & Lee 
2004). Obviously, in order to make such an analytical approach more useful for studying 
the evolution of subhalos and satellite galaxies, it is a prerequisite to use the correct initial 
phase space information of the infalling subhalos as found in the CDM cosmological model. 

Several authors have already studied the orbital parameters of the subhalos at the time 
of their mergers with the host halo (Tormen 1997; Vitvitska et al. 2002; Khochfar & Burkert 
2004; Benson 2005). Tormen (1997) and Khochfar & Burkert (2004) used high-resolution 
re-simulations of halos and identified the progenitors of these halos. The orbital distribution 
is then measured from the progenitors that are about to merge with the main progenitor. 
In contrast, Vitvitska et al. (2002) and Benson (2005) identified pairs of halos that are 
about to merge, and measured the orbital distribution of these pairs. Among these studies, 
Benson (2005) used a large set of cosmological simulations provided by the Virgo Consortium 
(Jenkins et al. 2001), and formed a large sample of such orbital pairs. From this, he presented 
fitting formulae for the distribution of the initial infall velocity of the subhalos at the virial 
radius of their host halo. 

Although the work of Benson (2005) is an important step forward to determine the initial 
condition of infall subhalos, there is an important question yet to be examined quantitatively, 
i.e., whether the mergers are isotropic in the position space or there is some preferential di- 
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rection for the mergers. Earlier studies suggest the accretion is anisotropic (e.g., Tormen 
1997). For example, Aubert, Pichon & Colombi 2004) found the infall takes place prefer- 
entially in the plane perpendicular to the direction defined by the spin of the halo. In this 
paper, we use a cosmological N-body simulation of 512 3 particles (Jing & Suto 2002; Kang 
et al. 2004) to investigate this problem in more detail. As we will show, the mergers are 
preferentially along the major axis of the host halos; we will quantify this anisotropic accre- 
tion as a function of the subhalo mass (see also Libeskind et al. 2005). Compared with the 
simulations used by Benson (2005), our simulation has higher force and mass resolutions, 
and the subhalos are better resolved. This makes it easier and more reliable to quantify the 
orbital distribution of subhalos. This is particularly important for small subhalos as they 
may lose their identities in low resolution N-body simulations when they are close to to the 
boundary of bigger halos. It is difficult to account for this population of missing subhalos if 
only one snapshot of the simulation output is used. 

The paper is structured as follows. In §2, we briefly describe the numerical simulations we 
use, and how the subhalos are identified. In §3, we study the physical properties of subhalos 
and present simple empirical fitting formulae to these properties. In §4 we summarize our 
main results and discuss areas for future improvement. 



2 N-BODY SIMULATION 

The simulation used in this paper is a P 3 M cosmological simulation of 512 3 particles in a 
box of lOO/i^Mpc. The cosmological model is the standard concordance model with the 
density parameter f2 mi0 = 0.3, the cosmological constant f^o = 0.7 and the Hubble constant 
h = Hq/(100 kms -1 Mpc -1 ) = 0.7. The initial density field is assumed to be Gaussian with a 
Harrison-Zel'dovich primordial power spectrum and with an amplitude specified by as = 0.9, 
where <t 8 is the r.m.s. fluctuation of the linearly evolved density field in a sphere of radius 
8 h^Mpc. This simulation, which started at redshift Zi = 72, is evolved by 5000 time steps 
to the present day (z = 0) with our vectorized parallel P 3 M code (Jing & Suto 2002). The 
force softening length r/f (S2 type, Hockney & Eastwood 1981) is 10 /i _1 kpc comoving, and 
the particle mass m p = 6.2 x 10 8 /i _1 M Q . Because these simulation parameters are very 
similar to those adopted in many high-resolution re-simulations of individual cluster halos 
(Moore et al. 1999, Jing & Suto 2000, Fukushige & Makino 2001, 2003, Power et al. 2003, 
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Diemand et al 2004), we have achieved a resolution that can resolve subhalos within massive 
host halos. 

The dark matter halos are identified at redshift z = in the simulation described above 
using the Friends-of-Friends method (FOF) with a linking length equal to 0.2 of the mean 
particle separation. The subhalos are then identified within the FOF halos with the SUBFIND 
routine (Springel et al. 2001). In Kang et al. (2004), the mass function of the subhalos 
was examined, and was found to be in good agreement with the subhalo mass functions 
obtained in previous halo re-simulations (Springel et al. 2001) down to a mass of about 
3.1 x 10 10 h~ 1 M & (50 particles). This indicates that the subhalos with more than iV su b = 50 
particles are resolved in the simulation. Considering the fact that subhalos can survive more 
easily in the outskirt than in the inner dense region of a host halo, we relax the lower limit 
to 20 particles for the subhalo mass. Table Q lists the number of particles we use to identify 
the subhalos and halos in several different combinations. And we only use the simulation 
output at z= to calculate the distribution. 

3 RESULTS 

3.1 The phase-space distribution in the spherical halo model 

We first consider the density profile of host halos as a sphere and search for subhalos within 
a spherical shell with radius r between 1 — Ar and 1 + Ar, where r is the distance between 
the centres of the subhalo and the host halo in units of the virial radius r vir of the host 
halo. We identified the "centres" of the halos as the lowest-potential particles using the 
SUBFIND routine (Springel et al. 2001; Kang et al. 2004); we also checked that the "centre" 
identified this way is very close to the centre of mass for most host halos. We normally take 
the thickness of the shell (Ar) to be 0.1 or 0.2 (see Table Q). The virial radius is determined 
according to the spherical collapse model (Kitayama & Suto 1996; Bryan & Norman 1999). 
In this section, we consider only those subhalos with an inward directed velocity. We treat the 
subhalos and host halos as point-mass particles, and determine the velocity of the subhalos 
at the time when their orbits first cross the virial radius of the host halo under gravity. We 
also compute the time t cr for each subhalo to cross from 1 + Ar to 1 — Ar according to their 
trajectories. Subhalos with a larger t cr will stay longer in the shell of 1 ± Ar than those 
with smaller t CT . During a time interval dt, the chance of observing this subhalo within the 
spherical shell will be dt/t cr . Thus we weight each subhalo by £~ x when we compute the 
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distributions of position, radial and tangential velocities. We will excise the subhalos that 
do not pass through one or both of the radial limits. 

The distributions of the radial (v r ) and tangential (vg) velocities of these subhalos are 
presented in Fig. 1; both velocities are in units of the circular velocity at the virial radius, 
f v i r . For this exercise, we only retain host halos with more than 500 particles and subhalos 
with more than 20 particles, corresponding to the selection parameter set A in Table d 
The v r distribution peaks around the virial velocity, while Vg peaks at a slightly smaller 
value, around 0.7. Both distributions drop essentially to zero beyond 1.5 virial velocity. 
The shape of these distributions qualitatively agrees with that obtained by Benson (2005). 
Quantitatively, however, both distributions appear slightly broader than those in Benson 
(2005). We do not know the exact reason for the difference, but it may be attributed to the 
different ways of selecting the subhalos. For low resolution simulations, subhalos close to the 
inner shell 1 — Ar may be difficult to identify with the friends-of-friends or the spherical 
overdensity methods. Benson (2005) attempted to correct for this effect by considering a 
minimum radius at which the subhalo can still be identified using his algorithm. With the 
SUBFIND routine (Springel et al. 2001), we are able to resolve the subhalos even when they 
are quite close to the host halos, so we do not need to make this complicated correction. 
Instead we simply consider all the subhalos within the spherical shell between radius 1 — Ar 
and 1 + Ar. 

A related quantity to v r and vg is the infall angle, defined as 

COS CKinfeu = -r^— n\T7? ■ K 1 ) 

( V r + v e) 1 

For a radially infalling subhalo onfall — 0. The distribution of o^aii is shown in the bottom 
left panel. The infall angle has a peak around 35°, and a full-width-at-half- maximum of 
about 50°. 

It is well known that halos accrete matter along the large-scale filaments that connect 
the halos, this implies that the merger of the subhalos into the host halos may be anisotropic 
(Lee, Jing, & Suto 2005). Since the halos are generally triaxial (Jing & Suto 2002), there 
could be a correlation between the direction of the subhalo merger and the shape of the 
host halo. We therefore determine the three principal axes for each halo from its inertial 
tensor within the virial radius. We define 9 to be the angle between the major axis of the 
host halo and the vector from the host halo centre to the centre of a subhalo, and 4> as the 
other polar angle (0 < < 2ir). Because of the limited sample size, we will focus on the 
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spatial distribution as a function of /i (= | cos#|), despite of the fact the distribution also 
depends on but more weakly than on 9. The probability distribution of fi, df/dn, for the 
subhalos is shown in the lower right panel of Fig. 1 . If the subhalos merge into the host halo 
isotropically, then we expect df /dfi to be unity. In contrast to this expectation, df/dfi of the 
subhalos increases strongly with fi, implying that subhalos are accreted more preferentially 
along the major axis of the host halos. 

Since dark halos are triaxial (Jing & Suto 2002), the dark matter density at the spherical 
virial radius is expected to be higher in the direction of the halo major axis. If the spatial 
distribution of subhalos follows that of the dark matter, we would expect a higher df / dfi in 
the direction of fi = 1. The question is whether the increase of df /dfi with /i can be fully 
explained by the shape of host halos. In the following subsection, we examine this question, 
i.e., consider the phase space distribution in the ellipsoidal coordinate system in the more 
realistic triaxial halo model. 



3.2 The phase-space distribution in the triaxial halo model 

For each host halo, we determine the axial ratios a/c and b/c from their inertia tensor within 
its virial radius r v ; r , where a, b, and c are the minor, median, and major axes of the halo 
respectively. If we rotate the coordinate (x,y,z) so that the new X, Y and Z coordinate 
axes are parallel to the minor, median, and major axes of the halo, the isodensity surfaces 
of a halo are approximately described by (Jing & Suto 2002) 
r 2 X 2 r 2 Y 2 

R 2 = -^r + ~ir + z 2 . (2) 



(i 



We define the boundary of an ellipsoidal halo R v [ r such that the total mass within R vir is equal 
to the virial mass of the halo in the spherical model. We then examine the orbital parameters 
and density of the subhalos at this surface. In analogy with the previous subsection, we take 
a shell of the upper and lower ellipsoidal radii (1 ± AR) in units of R V i T , and compute the 
distributions of the normal and tangential velocities (v n and v t ) relative to the ellipsoidal 
surface as well as the number density distribution. In this exercise, we only include the 
subhalos that have an inward directed velocity. The infall angle onfall can be similarly 
defined as in eq. ((TJ) with v r and vq replaced by v n and v t . 

The results are given in Fig. 2. The distributions of the normal and tangential velocities 
are very similar to those in spherical shells in the previous subsection. The subhalo number 
density n tr j(/x), which is the density averaged over <fi directions for a fixed \x on the surface 
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-R V ir, is shown as a function of /i in the middle right panel. Notice that in the triaxial halo 
model, a flat distribution in /i implies that the subhalos follow the idealised ellipsoidal shape. 
As expected, the dependence on \i becomes significantly weaker in the ellipsoidal shells, but 
there is still a preference for the subhalos to merge into their host halos along the major 
axis. To facilitate comparisons with previous studies, we also show in Fig. 2 (bottom two 
panels) the distributions of the total velocity v = (v n 2 + ft 2 ) 1 / 2 , and the orbital 'circularity' 
e = J / J c , where J is angular momentum and J c is the angular momentum of a circular orbit 
with the same energy. The orbital circularity has a broad peak around 0.5, indicating the 
infall is in neither purely radial nor tangential orbits. This result is in good agreement with 
previous studies (e.g., Tormen 1997, Fig. 4; Ghigna et al. 1998, Fig. 14). 

Next we examine the distribution of the velocities as a function of the angular positions of 
the subhalos relative to the major axis of the ellipsoidal host halo. To do this, we divide the 
ellipsoid into three equal bins in < // < 1. In Fig. 3, we show the distributions of the normal 
and tangential velocities for these three bins. The distributions are quite similar. The bottom 
two panels show the mean total velocity and the mean infall angle as a function of /i. Again 
these two quantities have little dependence on \i. Therefore we can use the distributions as 
shown in Figure 2 to more accurately describe the distribution of the subhalo velocities. 

We have also studied the distributions as a function of halo mass. For this purpose, 
we increase the lower limit of the host halos to 5000 particles (3.1 x 10 12 h~ 1 M Q ) which 
corresponds to the halo selection parameter set C in Table 1. The results are compared with 
those of the host halos with parameter set B, where we included host halos with more than 
500 particles. To increase the number of subhalos in the case of C, we have adopted Ar = 0.2. 
As shown in the two bottom panels of Figure 3, we found the statistical distributions change 
very little with reasonable changes of Ar. Figure 4 shows the results. No significant difference 
can be found for the two velocity distributions and the infall angle distribution. However, 
the number density dependence on /i becomes weaker when we increase the lower limit of 
the host halo mass and hence include more subhalos with smaller M su b/Mh os t- If the smaller 
subhalos follow more closely the triaxial density profile, then the trend seen in Fig. 4 can be 
understood because the mass density should follow the triaxial density model. 

To see the preferential accretion of more massive subhalos along the major axis, in Fig. 
5 we divide the subhalos in six bins of M su b/Mh os t (with roughly equal numbers) and show 
their distributions of \x. The mass bins are listed in Table 3. The lowest mass-bin subhalos 
have only a weak dependence on /i, indicating that they are roughly consistent with the 
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triaxial shape. However, the highest mass subhalos (top left panel) with M su b/Mh os t ~ 0.08 
shows dramatic deviation from from the ellipsoidal accretion, with a strong peak around 
/i = 1 (i.e., 9 = 0°). Fig. 5 clearly illustrates that more massive subhalos are accreted with 
a much stronger preference along the major axis of the host halo. This may have observable 
consequences for the satellite galaxies we see today; we return to this important point briefly 
in the discussion. 

We have shown that the infall of subhalos are preferentially along the major axis of the 
parent halo determined by the moment of inertia. An alternative way of defining the orien- 
tation of a halo is using its spin axis, which can be defined by the total angular momentum 
of the particles within the virial radius. Aubert, Pichon & Colombi (2004) found that the 
infall takes place preferentially in the plane perpendicular to the direction defined by the 
spin of the halo. Fig. El shows the distribution of the angle between the spin axis and the line 
connecting the centres of the parent halo and the satellite in our simulation. Clearly there 
is a (weak) preference for satellites to be perpendicular to the spin axis, in agreement with 
the findings of Aubert et al. (2004, see their Fig. 8, which is based on a more sophisticated 
spherical harmonics analysis). This is to be expected as the spin axis is statistically perpen- 
dicular to the major axis of the halo, as shown by previous studies (e.g. Warren et al. 1992, 
Dubinski 1992). 

Many analytical models of the subhalo population require an accurate knowledge of the 
initial conditions for the subhalos. In the following, we provide empirical fitting formulae 
found in our numerical simulations. To facilitate comparisons with Benson (2005), we adopt 
the same functional forms to fit the velocity distribution, specifically, 

f(v n ,v t ) = aiv t exp[-a 2 (v t - a 9 ) 2 - &i(-u t ){v n - b 2 (v t )} 2 ] (3) 
where 

bi(v t ) = a 3 exp[-a 4 (ft - a 5 ) 2 ], (4) 

b 2 {v t ) = a 6 exp[-a 7 (f t - a 8 ) 2 ]. (5) 

The fitted curves are shown in Fig. 2, and the fit parameters are listed in Table El We also 
use the following function 

df /dfi = p + pi exp(p 2 Z^ 2 ), = cos 6 1 , < /i < 1 (6) 

to describe the angular distribution of the subhalo accretion. The function has two param- 
eters, pi and p 2 , and po is chosen such that the function is properly normalised when \i is 
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integrated from to unity. The fit parameters are listed in Table As can be seen, the 
fitting functions match our simulation results quite well. 

4 DISCUSSIONS 

We have used a high-resolution simulation to study the initial conditions of subhalos when 
they merge into their host halos. Most of our results are in good agreement with Benson 
(2005). One finding of our study is that massive subhalos are accreted more preferentially 
along the major axis of the host halos than the less massive ones. Our subhalos are identified 
at present day [z = 0), but the same trend for massive subhalos should apply at high 
redshifts (Kravtsov et al. 2004). If the more massive subhalos house satellite galaxies and 
they survive until the present day, then the satellites should show a more planar geometry 
along the major axis. Interestingly, the satellite galaxies in the Milky Way appear to lie in 
a great disk (Kroupa et al. 2005; see also Willman et al. 2004), almost perpendicular to the 
stellar disk. If the major axis of the dark matter halo is perpendicular to the stellar disk in 
the Milky Way, then such a distribution, while puzzling at first, is naturally expected in the 
CDM (Kang et al. 2005; Libeskind et al. 2005; Zentner et al. 2005b). 

Our result that more massive subhalos are more preferentially accreted along the major 
axis of the host halo can be understood in the cosmic web theory (Bond, Kofman, & Pogosyan 
1996; Lee, Jing, & Suto 2005). In this theory, the coherence of the initial tidal field forms 
one dimensional filaments. The halos are bridged by the filaments, and the merger of halos 
is expected to occur along the filaments. If the major axis of host halos is determined by the 
neighboring massive filaments, we would expect subhalos to merge preferentially along the 
major axis. A quantitative computation is possible for the distribution of df/dfi following 
Lee et al. (2005) which will be explored in a future work. 

The initial conditions we found can be used in analytical formalisms (e.g., Oguri & Lee 
2004; Taylor & Babul 2004) to predict the subsequent evolutions of subhalos and compare 
these with observations. Our results make it straightforward to generate Monte Carlo real- 
izations of the subhalo population. The extended Press & Schechter formalism accurately 
predict the mass function of the subhalos. For a given halo, the triaxial halo shape can be 
sampled using Monte Carlo with the fitting formulae given by Jing & Suto (2002). The 
velocities of the halo can then be obtained using eq. (jSJ). To describe the angular position 
of the satellite, the polar angles (9, 0) are needed. To the first order, the polar angle (4>) is 
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approximately uniformly distributed between to 2tt. The angle 9 can be generated using 
eq. JHJ). The mass dependence of the angle 9 can be properly taken into account by choosing 
different fit parameters (in Table 3) depending on the value of M su b/Mh os t- 

In this paper we have analyzed the phase-space information of subhalos located in a thin 
shell close to the virial radius. In practice, we found that about 30% of subhalos within the 
shell have outgoing velocities. Some of these subhalos and perhaps even some inward moving 
subhalos could have been inside the host halo already, thus some of these subhalos could 
have evolved dynamically within the host halo and their phase-space information cannot 
be regarded as "initial" conditions. The number of secondary infall halos can be crudely 
estimated as follows. Let us approximate the parent halo as spherical, and assume that 
there are N subhalos falling into the parent halo for the first time. Due to symmetry, we 
can assume all subhalos move along the positive x-axis, and a fraction / of the subhalos 
will emerge as outgoing subhalos on the opposite side of the halo. Presumably all these 
fN subhalos will be accreted again in a second infall, and some fraction (presumably £ f) 
of these subhalos will emerge as outgoing ones on the positive x side. So in total there are 
N + fN infall subhalos, and fN+f' 2 N outgoing subhalos. The ratio of infalling and outgoing 
subhalos is therefore (N + f'N) /(f'N + f' 2 N) pa 7/3, implying / ~ 3/7. Hence the fraction of 
the secondary infall subhalos is about fN/ (N + fN) £ 30%. Gill et al. (2005) investigated 
satellite galaxies in the outerskirts of galaxy clusters from their high-resolution simulations 
(see also Gill et al. 2004a,b). They found that approximately one half of the galaxies with 
current cluster-centric distance in the interval 1-2 virial radii of the host are 'backsplash' 
galaxies that once penetrated the cluster potential. If one half of these 'backsplash' galaxies 
have negative radial velocities, then this suggests that the fraction of subhaloes undergoing 
secondary infall is about ~ 25%, consistent with our estimate above. 

To more properly account for this population of evolved subhalos in detail, one must 
track the subhalos as a function of redshift. This will also be useful for studying the time 
evolution of initial conditions. Another issue not studied in detail here is the <fi dependence 
of the subhalo spatial distribution. This quantity is taken to be uniform between to 27r, 
while approximately correct to the first order, a more realistic treatment is desirable. We 
plan to return to these issues in future works. 



12 Wang, Jing, Mao & Kang 
ACKNOWLEDGMENTS 

We thank the referee for a helpful report which improved the paper. The research is sup- 
ported by NKBRSF (G19990754), NSFC (Nos. 10125314, 10373012), and Shanghai Key 
Projects in Basic research (No. 04jcl4079). SM acknowledges the financial support of Chi- 
nese Academy of Sciences and the European Community's Sixth Framework Marie Curie 
Research Training Network Programme, Contract No. MRTN-CT-2004-505183 "ANGLES". 
HYW and SM wish to thank the hospitality of Shanghai Astronomical Observatory during 
several scientific visits. 



The phase-space distribution of inf ailing dark matter subhalos 13 

REFERENCES 

Aubert D. Pichon C. Colombi S., 2004, MNRAS, 352, 376 
Benson A. J., 2005, MNRAS, 358, 551 

Benson A. J., Lacey C. G., Frenk C. S., Baugh C. M., Cole, S. MNRAS, 2002, 351, 1215 

Bertone C, Hooper D., Silk J., 2004, Phys. Rep. 405, 279 

Bond J., R., Kofman L., Pogosyan D., 1996, Nature, 380, 603 

Bryan G. L., Norman M. L., 1998, ApJ, 495, 80 

Colin P, Klypin A. A., Kravtsov A. V., 2000, ApJ, 539, 561 

Diemand J., Moore B., Stadel J. 2004, MNRAS, 352, 535 

Dubinski J., 1992, ApJ, 401, 441 

Gao L., De Lucia G., White S. D. M., Jenkins A., 2004a, MNRAS, 352, LI 

Gao L., White S. D. M., Jenkins A., Stoehr F., Springel V., 2004b, MNRAS, 355, 819 

Ghigna S., Moore B., Governato F., Lake G., Quinn T., Stadel J., 1998, MNRAS, 300, 146 

Ghigna S., Moore B., Governato F., Lake G., Quinn T., Stadel J., 2000, ApJ, 544, 616 

Gill S. P. D., Knebe A., Gibson B. K. 2004a, MNRAS, 351, 399 

Gill S. P. D., Knebe A., Gibson B. K., Dopita A. D., 2004b, MNRAS, 351, 410 

Gill S. P. D., Knebe A., Gibson B. K., 2005, MNRAS, 356, 1327 

Hockney R. W., Eastwood J. W., 1981, Computer Simulation Using Particles (New York: 
McGraw-Hill) 

Jing Y. P., Suto Y., 2002, ApJ, 574, 538 

Kang X., Jing Y. P., Mo H. J., Bonier G., 2004, ApJ, in press ( |astro-ph/0408475D 

Kang X., Mao S., Gao L., Jing Y. P., 2005, A&A, 437, 383 

Kitayama T., Suto Y., 1996, ApJ, 469, 480 

Klypin A., Kravtsov A. V., Valenzuela O., 1999, ApJ, 522, 82 

Kravtsov A. V., Gnedin O. Y., Klypin A. A., 2004, ApJ, 609, 482 



Khochfar S., Burkert A., 2004, MNRAS, submitted jastro-p h/0309611) 
Kroupa P., Theis C, Boily C. M., 2005, A&A, 431, 517 
Kochanek C.S., Dalai N., 2004, ApJ, 610, 69 
Lee J., Jing Y.P., Suto Y., 2005, ApJ, in press 

Libeskind N. L, Frenk C. S., Cole S., Helly J. C, Jenkins A., Navarro J. F., Power C, 
2005, |astro-ph/0503400 



Mao S., 2004, IAU 220, Eds. M. A. Walker, & K. C. Freeman. ASP 237, p85 



14 Wang, Jing, Mao & Kang 

Mao S., Schneider P., 1998, MNRAS, 295, 587 

Mao S., Jing Y. P., Ostriker J. P., Weller J. 2004, ApJ, 604, L5 

Moore B., Ghigna S., Governato F. et al, 1999, ApJ, 524, L19 

Oguri M., Lee J., 2004, MNRAS, 355, 120 

Okamoto T., Habe A., 1999, ApJ, 516, 591 

Sheth R. K., 2003, MNRAS, 345, 1200 

Springel V., White S. D. M., Tormen G., Kauffmann G., 2001, MNRAS, 328, 726 
Spergel D.N.S. et al., 2003, ApJS, 148, 175 

Stoehr F., White S. D. M., Springel V., Tormen G., Yoshida N. 2003, MNRAS, 345, 1313 
Taylor J. E., Babul A., 2001, ApJ, 559, 716 
Taylor J. E., Babul A., 2004, MNRAS, 348, 811 
Taylor J. E., Silk J., 2003, MNRAS, 339, 505 
Tormen G., 1997, MNRAS, 290, 411 

Vitvitska M., Anatoly A., Kravtsov A. V., Wechsler R. H., Primack J. R., Bullock J. S., 
2002, ApJ, 581, 799 

Warren M. S., Quinn P. J., Salmon J. K., Zurek W. H., 1992, ApJ, 399, 405 

Willman B., Governato F., Dalcanton J. J., Reed D., Quinn T. 2004, MNRAS, 353, 639 

van den Bosch F. C, Lewis G. F., Lake G., Stadel J., 1999, 515, 50 

van den Bosch F. C, Tormen G., Giocoli C, 2005, MNRAS, 359, 1029 

White S. D. M., Rees M. J., 1978, MNRAS, 183, 341 

Zentner A. R., Bullock J. S., 2003, 598, 49 

Zentner A. R., Berlind A. A., Bullock J. S., Kravtsov A. V., Wechsler, R. H., 2005a, ApJ, 
624, 505 

Zentner A. R., Kravtsov A. V., Gnedin O. Y., Klypin A. A. 2005b, astro-ph/0502496 



The phase-space distribution of infalling dark matter subhalos 15 



Table 1. Halo and subhalo selection parameters: iVhalo an d N S ub are the lower limits of the particle number for the host halos 
and subhalos respectively. The subhalos are selected within a spherical shell from 1 — Ar to 1 + Ar in units of the virial radius. 
The numbers of selected subhalos are listed in the last row. 



Model 


A 


B 


C 




500 


500 


5000 




20 


20 


20 


Ar 


0.1 


0.2 


0.2 


number of subhalos 


2606 


4437 


3470 



Table 2. Fit parameters for the velocity distributions as defined in cq. 13*1 , 



parameter A B 



01 






0.1 


2.12 


2.62 


0.3 


2.90 


4.48 


(14 


-0.333 


-0.525 


«5 


-0.490 


0.238 


Og 


1.17 


1.20 


07 


0.155 


0.140 


a* 


-0.564 


-0.731 


ag 


0.314 


0.294 



Table 3. Fitting parameters for the density distribution n tI i(/i) in different mass bins as defined in eq. JSJ. Each bin is specified 
by the lower M\ and upper M u limits of the subhalo mass in units of the mass of the host halo. 





all subhalos 


bin 1 


bin 2 


bin 3 


bin 4 


bin 5 


bin 6 


M u 


1 


1 


0.07943 


0.03946 


0.02542 


0.01709 


0.00834 


Mi 


0.0 


0.07943 


0.03946 


0.02542 


0.01709 


0.00834 


0.0 


pO 


0.644 


0.289 


0.259 


0.361 


0.377 


-4.21 


-617.8 


pl 


0.194 


6.33X10" 4 


7.37X10" 2 


0.159 


0.258 


4.98 


618.8 


p2 


1.494 


9.953 


4.281 


2.919 


2.040 


0.130 


2.505 xl0~ 4 
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Figure 1. Orbital parameters and spatial distributions for the subhalos selected according to the selection parameter set A 
in Table H\ in spherical halo model. The upper left and right panels show the distributions of radial and tangential velocities, 
respectively. The lower left panel shows the distribution of the infall angle, denned in eq. 0. In the lower right-hand panel, 
fill square shows the distribution of the number density of subhalos per unit area as a function of the angle 9 between the 
vector from the subhalo centre to the host halo centre and the major axis of the host halo. For an isotropic population, the 
distribution df/dfj, is a horizontal line with amplitude unity. We also show the cumulative distributions for the simulation and 
an isotropic distribution (open triangles and open circles) in the lower right panel. 
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Figure 2. Orbital parameters and spatial distributions for the subhalos selected according to the parameter selection set A in 
Table H\ in triaxial halo model. The upper left and right panels show the distributions of the normal and tangential velocities 
respectively. The middle left panel shows the distribution of the infall angle. The middle right panel shows the number density 
n tI i(fi) of subhalos at the surface i? v j r as a function of 9, where 8 is the angle between the major axis of the host halo and the 
vector from the subhalo centre to the host halo centre. The fitting curves as given in eqs. J3J and JHJ are shown as a thin line 
in each panel. The bottom left and right panels show the orbital circularity and the total infall velocity respectively. Notice 
that the vertical scale in the middle right panel is different from that in the bottom right panel in Fig. 1. 
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Figure 3. The upper panels plot the distributions of the normal and tangential velocities in different ranges of fi = cos0, 
where 6 is the angle between the vector from the subhalo centre to the host halos centre and the major axis of the host halo. 
The lower left panel shows the average velocity while the lower right panel shows average infall angle as a function of fj.. 
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Figure 4. The upper panels show the distributions of the normal and tangential velocities respectively in the triaxial halo 
model for different subhalo masses (Sets B and C). The lower left panel shows the distribution of the infall angle (ai n f a ii) while 
the lower right panel shows the number density n tr ;(/i) of subhalos at the surface -R v i r , where fi = cos 6 and 8 is the angle 
between the major axis of the host halo and the vector connecting the centres of the subhalo and the host halo. The difference 
between the models B and C in the bottom right panel is due to the different lower limits of the host halo mass (see §3.2). 
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Figure 5. The number density n tr i(p-) of subhalos at the surface i? v ir together with the best-fit curves in the triaxial halo 
model. The subhalos are divided into six bins of M su ]j/Mh ost , and the result for each bin is presented in one panel. The mass 
range for each bin and the corresponding fit parameters (as defined in eq. |5) are listed in Table 151 The distribution for all the 
subhalos is shown in the top left panel. Notice that the vertical scales are different in different panels. 
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Figure 6. The distribution of fi spln = cos 6 sp ; n in the tri-axial model, where S pin is the angle between the spin axis and the 
line connecting the centres of the parent and satellite halos. The straight line shows the prediction for an isotropic distribution. 



